o 

(N 






Ph. 



O 



o 
o 



Acta Mechanica Sinica manuscript No. 

(will be inserted by the editor) 



Tianshu Luo Yimu Guo 

An Examination of the Time- Centered 
Difference Scheme for Dissipative 
Mechanical Systems from a Hamiltonian 
Perspective 

Received: date / Accepted: date 



Abstract In this paper, we have proposed an approach to observe the time-centered 
difference scheme for dissipative mechanical systems from a Hamiltonian per- 
spective and to introduce the idea of symplectic algorithm to dissipative systems. 
The dissipative mechanical systems discussed in this paper are finite dimensional. 
This approach is based upon a proposition: for any nonconservative classical me- 
chanical system and any initial condition, there exists a conservative one; the two 
systems share one and only one common phase curve; the Hamiltonian of the 

rS^ • conservative system is the sum of the total energy of the nonconservative sys- 

j^ I tem on the aforementioned phase curve and a constant depending on the initial 

condition. Hence, this approach entails substituting an infinite number of conser- 
vative systems for a dissipative mechanical system corresponding to varied initial 
conditions. One key way we use to demonstrate these viewpoints is that by the 
Newton-Laplace principle the nonconservative force can be reasonably assumed 
^1 to be equal to a function of a component of generalized coordinates qj along a 

phase curve, such that a nonconservative system can be reformulated as countless 
conservative systems. The advantage of this approach is such that there is no need 

»vj i to change the definition of canonical momentum and the motion is identical to 

that of the original system. Therefore, first we utilize the time-centered difference 
scheme directly to solve the original system, after which we substitute the numer- 
ical solution for the analytical solution to construct a conservative force equal to 
the dissipative force on the phase curve, such that we would obtain a substitut- 
ing conservative system numerically. Finally, we use the time-centered scheme 
to integrate the substituting system numerically. We will find an interesting fact 
that the latter solution resulting from the substituting system is equivalent to that 
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of the former. Indeed, there are two transition matrices within time grid points: 
the first one is unsymplectic and the second symplectic. In fact, the time-centered 
scheme for dissipative systems can be thought of as an algorithm that preserves 
the symplectic structure of the substituting conservative systems. In addition, via 
numerical examples we find that the time-centered scheme preserves the total en- 
ergy of dissipative systems. According to such behaviors, we might explain why 
some algorithms, e.g., the time-centered Euler scheme, are better than other un- 
symplectic algorithms for dissipative systems, such that we might choose better 
algorithms and introduce the idea of this paper to more symplectic algorithms. 

Keywords Hamiltonian, dissipation, non-conservative system, damping, 
symplectic algorithm 



1 Introduction 

Feng ll3]fT4lfT5l l4l . Zhong |fT7]l 6l[T3l and Marsden|[9| have developed a series of sym- 
plectic algorithms for conservative systems and proposed common theories for the 
construction of symplectic algorithms. 

Feng||3l has investigated some existing old numerical algorithms from a Hamil- 
tonian perspective, such as the simplest symplectic algorithm, the Euler time- 
centered difference scheme. He have explained why for a conservative system the 
Euler time-centered difference scheme is more accurate than other unsymplectic 
schemes from a Hamiltonian perspective. Because the time-centered difference 
scheme is derived from Hamilton's equation, the algorithm can preserve symplec- 
tic structure and mechanical energy. The mechanical energy-preserving behavior 
was proven by Xing ||I6l . who found that for dissipative problems this algorithm 
has good mechanical energy-preserving characteristics. We will state the reason 
using the Hamiltonian description of dissipative mechanical systems, which are 
finite dimensional in this paper, and we will then apply the idea of symplectic 
algorithms to the analyses of the time-centered difference scheme for dissipative 
systems. Nevertheless, since Hamilton originated Hamilton equations of motion 
and Hamiltonian formalism, it has been stated in most classical textbooks that the 
Hamiltonian formalism focuses on solving conservative problems. 

If one needs to apply symplectic algorithms to dissipative systems, one must 
convert a dissipative system into a Hamiltonian system or find some relationship 
between the dissipative system and a conservative one. An attempt[7| has been 
made to apply symplectic algorithms directly to dissipative systems. The transi- 
tion matrix between phase variable vectors must be unsymplectic. Symplectic al- 
gorithms called variational integrators Q were derived from discretization of the 
variational principle for conservative systems 

Liq{t),q{t))dt = 0, 

which is equivalent to Hamilton's equation. Therefore, the algorithms called varia- 
tional integrators are natively symplectic. But the direct variational integrators for 
dissipative systems were derived from discretization of the Lagrange-d' Alembert 
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principle 

rb rb 

L{q{t),q{t))At+ / F{q{t),q{t))8qAt=0, 
Ja 

where F{q(t),q(t)) is a nonconservative force. Obviously Hamilton's equation 
cannot be derived from the Lagrange-d' Alembert principle. ZhongQ considered 
that discretization of Hamihon's equations or the Hamihonian least action vari- 
ational principle can lead to symplectic transition matrices. Therefore, the varia- 
tional integrators applied directly to dissipative systems cannot be considered as 
symplectic schemes. 

ZhongtlTl attempted to convert a damping Duffing equation into a conserva- 
tive system by redefining the position variable q and the canonical momentum p 
and then applying the time-fem method to this conservative system. His method 
is to multiply q by the reciprocal of the amplitude decay coefficient such that the 
momentum is redefined. The characteristic of the original dissipative system has 
changed entirely. 

Although several other approaches have been proposed to represent dissipative 
systems as Hamiltonian formalism, these approaches might not be accepted by re- 
searchers in geometrical mechanics. For instance, Morrison lllOllllI and Salmon lll2ll 
focused on the conservative system or some special dissipative systems, e.g. an os- 
cillator with gyroscopic damping. Morrison 1 101 wrote, 'The ideal fluid description 
is one in which viscosity or other phenomenological terms are neglected. Thus, as 
is the case for systems governed by Newton's second law without dissipation, 
such fluid descriptions posses Lagrangian and Hamiltonian descriptions.' If there 
had existed an approach appropriate for representing an oscillator with damping 
as Hamiltonian formalism, these researchers would have attempted to extend the 
Hamiltonian description to other dissipative problems. 

Marsden [8| and other researchers applied the equations as below to the prob- 
lem of stability of dissipative systems 

dH pfdr 
dqt \dqi 

i^ = §^, (1) 

where the position vector r depends on the canonical variable {q^p}, i.e. r{q,p), H 
denotes Hamiltonian, and F{dr/dqi) denotes a generalized force in the direction 
i, i = 1, . . . ,n. Marsden considered that Eqs.([Tll was composed of a conservative 
part and a non-conservative part. Eq.([T| apparently is not a Hamilton's equation 
but only a representation of dissipative mechanical systems in the phase space. Al- 
though one can utilize the approaches discussed in some papers to convert Eq.([T]) 
into a Hamiltonian system, one must first change the definition of the canonical 
momentum of the system. If one uses symplectic algorithms to solve the Hamil- 
tonian system, the numerical solution will lose the physical characteristics of the 
original system, because the phase flow of the original system is different from 
that of the new system. We need a Hamiltonian system that shares common phase 
flow or solution with the original system. But this demand cannot be satisfied, be- 
cause it conflicts with Louisville's theorem. Therefore, we would have to attempt 
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to find some relationship between dissipative systems and conservative ones, such 
that we can introduce the concept of symplectic algorithms to dissipative problems 
and explain the time-centered difference scheme from Hamiltonian viewpoints. 

Based on Eq.([T]), in this paper we will attempt to demonstrate that a dissipa- 
tive mechanical system shares a single common phase curve with a conservative 
system. In the light of this property, we will propose an approach to substitute a 
group of conservative systems for a dissipative mechanical system. In the follow- 
ing section, we will illustrate the relationship between a dissipative mechanical 
system and a conservative one. 



2 Relationship between a Dissipative Mechanical System and a Conservative 
One 

2.1 A Proposition 

Under general circumstances, the force F is a damping force that depends on the 
variable set qi,--- ,q„,qi,--- ,q„. Fi denotes the components of the generalized 
force F. 

Fiiqi,--- ,qn,qi,--- ,qn)=F i-^j ■ (2) 

Thus we can reformulate Eq.([Tll as follows: 

Pi = -3— +AW1,--- ,q,uq\,--- An) 
oqi 

q. = §^. (3) 

Suppose the Hamiltonian quantity of a conservative system without damping is 
H. Thus we may write a Hamilton's equation of the conservative system : 

dH 
dqi 

1i = ^-- (4) 

aPi 

We do not intend to change the definition of momentum in classical mechanics, 
but we do require that a special solution of Eq.© is the same as that of Eq.OJ. 
We may therefore assume a phase curve /of Eq.© coincides with that of Eq.(|4]l. 
The phase curve y corresponds to an initial condition qio,Pio- Consequently by 
comparing Eq.OJ and Eq.©, we have 

- Fi{qi,--- ,q„,qi,--- ,qn)\y, 
r 

(5) 



dH 
dqt 


dH 

y dqt 


dH 


dH 


dpi 


y dPi 
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where 4^ 



dH 



did 



and 



dH 



denote the values of these partial derivatives 
, denotes the value of the force 



on the phase curve /and Fi{qi,--- ,qn,qi,-- ■ ■,qn)\y' 
Fj on the phase curve y. In classical mechanics the Hamiltonian // of a conserva- 
tive mechanical system is mechanical energy and can be written as: 



H = 



\dqij ' jy\dpi) 



const I , 



(6) 



where const\ is a constant that depends on the initial condition described above. 
If qi = 0,pi = 0, then const\ = 0. The Einstein summation convention has been 
used this section. Thus an attempt has been made to find 77 1 through line integral 
along the phase curve y of the dissipative system 

AnA\-:--- An) 



Jy dqi ' Jy 



dH 



dqi 



y oPi Jy opi 

Analogous to Eq.©, we have 



H\ 



dH_ 

dqi 



dqi 



dH 
dp 



■dpi + const2, 



(7) 



(8) 



where const2 is a constant which depends on the initial condition. Substituting 
Eq.(|6]l(l7]l into Eq.®, we have 



H\ 



H / Fi{qi,--- ,q„,qu--- ,q,i)<iqi + const. 



(9) 



where const = const2 — const i, and H = H\y because H is mechanical energy 
of the nonconservative systemO]). According to the physical meaning of Hamilto- 
nian, const \, const2 and const are added into Eq.©®® respectively such that the 
integral constant vanishes in the Hamiltonian quantity. Arnold[|2| had presented 
the Newton-Laplace principle of determinacy as, 'This principle asserts that the 
state of a mechanical system at any fixed moment of time uniquely determines all 
of its (future and past) motion.' In other words, in the phase space the position 
variable and the velocity variable are determined only by the time f. Therefore, 
we can assume that we have already a solution of Eq.([3]l 



1i 

4i 



qi{t) 
4i{t), 



(10) 



where the solution satisfies the initial condition. We can divide the whole time 
domain into a group of sufficiently small domains and in these domains qi is 
monotone, and hence we can assume an inverse function t = t{qi). If f = f(^/) 
is substituted into the nonconservative force F,L, we can assume that: 



Fi{qi{t{qi)),- ■ ■ ,qn{t{qi))A\ {t(qi)), ■ ■ ■ An{t{qi)))L = ^i{qi) 



(11) 
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where ^,- is a function of qi alone. In Eg. (II lb the function Fj is restricted on the 
curve y, such that a new function ^i{qi) yields. Thus we have 

JFiAqt = r ^i{qi)Aqi = Wi{qi) - Wi{qio). (12) 

■'7 Mia 

According to Eq. (ll2b the function ^, is path independent, and therefore J?, can be 
regarded as a conservative force. For that Eg. dill) represents an identity map from 
the nonconservative force F on the curve y to the conservative force J^, which is 
distinct from Ft. Eq. dllb is tenable only on the phase curve y. Consequently the 
function form of ^, depends on the aforementioned initial condition; from other 
initial conditions ^, with different function forms will yield. 

According to the physical meaning of Hamiltonian, const is added to Eq.(|9]l 
such that the integral constant vanishes in Hamiltonian quantity. Hence const — 
~W,(^,o). Substituting Eq. (ll2b and const = — W,(^,o) into Eq.(|9l), we have 

H\^ = H-Wi{qi) (13) 

where — W,(^,) denotes the potential of the conservative force ^, and W/(g,) is 
equal to the sum of the work done by the nonconservative force F and const. In 
Eq. dOl) H and H are both functions of qi and W,(^,) a function of ^,. Eq. dl3l) 
and Eq.(|9]l can be thought of as a map from the total energy of the dissipative 
system© to the Hamiltonian of the conservative system®. Indeed, 77 1 and the 

total energy differ in the constant const = —Wi{qio). When the conservative system 
takes a different initial condition, if one does not change the function form of H I , 

one can consider H | as a Hamiltonian quantity H, 

H = H\^ = H~Wi(qi) (14) 

and the conservative system© can be thought of as an entirely new conserva- 
tive system. 

Based on the above, the following proposition is made: 

Proposition 1 For any nonconservative classical mechanical system and any ini- 
tial condition, there exists a conservative one; the two systems share one and only 
one common phase curve; the value of the Hamiltonian of the conservative sys- 
tem is equal to the sum of the total energy of the nonconservative system on the 
aforementioned phase curve and a constant depending on the initial condition. 

Proof First we must prove the first part of the Proposition[Tl i.e. that a conservative 
system with Hamiltonian presented by Eq. dl4l) shares a common phase curve with 
the nonconservative system represented by Eq.©. In other words the Hamiltonian 
quantity presented by Eq. (ll4b satisfies Eq.(|5]l under the same initial condition. 
Substituting Eq. (ll4b into the left side of Eq.(|5]l, we have 

dH{qi,Pi) _ dHiqi,pi) dWjiqj) 

(15) 



dqi dqi dqi 

dHjqi^Pi) ^ dH{qi,Pi) dWjjqj) 

dpi dpi dpi 
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It must be noted that although qi and /?, are considered as distinct variables in 
Hamilton's mechanics, we can consider q, and qi as dependent variables in the 
process of constructing ofH. At the trajectory 7 we have 

dqi dqi 

dWMj) 



dpi 



0, (16) 



where ^({qi) is equal to the damping force Ft on the phase curve y. Hence un- 
der the initial condition qo,po, Eq.(|5]l is satisfied. As a result, we can state that 
the phase curve of Eq.(|4]l coincides with that of Eq.© under the initial condi- 
tion; and H represented by Eq. dHI) is the Hamiltonian of the conservative system 
represented by Eq.Q. 

Then we must prove the second part of Proposition [T] the uniqueness of the 
common phase curve. 

We assume that eq.Q shares two common phase curves, 71 and 72, with eq.©. 
Let a point of Y\ at the time f be zi , a point of 72 at the time t Z2, and g' the Hamil- 
tonian phase flow of eq.(|4|. Suppose a domain fl at t which contains only points 
zi and Z2, and £2 is not only a subset of the phase space of the nonconservative 
system© but also that of the phase space of the conservative system©. Hence 
there exists a phase flow g' composed of Y\ and 72, and g' is the phase flow of 
eq.© restricted by Q. According to the following Louisville's theoremUl: 

Theorem 1 The phase flow of Hamilton 's equations preserves volume: for any 
region D we have 

volume of g'D = volume of D 

where g' is the one-parameter group of transformations of phase space 

g':{p{0),qm^:ip(.t),qit)) 



g' preserves the volume of £2. This implies that the phase flow of eq.© g' pre- 
serves the volume of Q too. But the system © is not conservative, which conflicts 
with Louisville's theorem; hence only a phase curve of eq.© coincides with that 
of eq.©. 

n 



2.2 An Example in Vibration Mechanics 

Take an «-dimensional oscillator with damping as an example, the governing equa- 
tion of which is as below: 

q+Cq+Kq = 0, (17) 
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where q=[qi,...,q„] , superscript T denotes a matrix transpose, 



C = 



Cii . 


■ ■ Cin 


,K = 


Cnl ■ 


• • ^nn 





K\\ ... K 



12 



K21 . . . K22 



, and Cij and Kjj are constants. 

It is complicated to solve Eq.( ll7l) . If Eq. JTTT l is higher dimensional, it is almost 
impossible to solve Eg. ( 117b analytically. Therefore we assume that a solution ex- 
ists already. 



q = q{t) = [qi{t),...,q„{t)]. 
Suppose a group of inverse functions 

t = t{qi),...,t = t{qn). 



(18) 



(19) 



As in Eq. dlll) we can consider that the damping forces are equal to some conser- 
vative force under an initial condition 



cii^l =Pll(^l) ■•• ci„q„ = pi„{qi) 



Cniqi = P2\{qn) ■ ■ ■ Cnnqn = Pnn{qn) ■ 



(20) 



For convenience, these conservative forces can be thought of as elastic restoring 
forces: 



Pii(^i) = k:ii(^i)^i ... p\n{q\) = K\n{q\)q\ 

Pn\{q\) = Ki\{qn)qn ■■■ Pnniq,,) = Knn{.qn)qn- 

An equivalent stiffness matrix K is obtained, which is a diagonal matrix 

n 

Ku = Y^Kn{qi). 



(21) 



/=i 



Consequently an n-dimensional conservative system is obtained 
q+{K+k)q = 



(22) 



(23) 



which shares a common phase curve with the «-dimensional damping svstem dTTl l. 
The Hamiltonian of Eqs.( l23l) is 



H = ^q^q + ^q''Kq+ l\Kqfdq, 



(24) 



where is a zero vector. H in Eq.( l24b is the mechanical energy of the conserva- 
tive svstem(l23]l, because jQ{Kq)^dq is a potential function such that^ doest not 
depend on any path. 
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A Dissipative Meclianical System with Varied Initial Conditions 



Conservative System 1 



Conservative System 2 



Conservative System n 



Fig. 1 A Dissipative Mechanical System and Conservative Systems 



2.3 Discussion 



Based on the above, we can outline the relationship between a dissipative me- 
chanical system and a group of conservative systems by means of Fig. [T] The 
relationship can be stated from two perspectives: 

If one explains the relationship from a geometrical perspective, one can ob- 
tain Proposition [T] In this paper the conservative systems (|4]l and (123 b are called 
the substituting systems. Although a substituting system shares a common phase 
curve with the original system, under other initial conditions the substituting sys- 
tem exhibits different phase curves. Therefore the phase flow of the substituting 
system differs from that of the original system, it follows that the substituting sys- 
tems is not equal to the original system. According to Louisville's theorem ([T|, the 
phase flow of the original dissipative system Eq.lO certainly does not preserve its 
phase volume, but the phase flow of the substituting conservative Eq.(|4]l does. 

One also could explain the relationship from a mechanical perspective. It is 
known that there are non-conservative forces in a nonconservative system. The 
total energy of the nonconservative system consists of the work done by noncon- 
servative forces. Hence the function form of the total energy depends on a phase 
curve i.e. under an initial condition. If one constrains the total energy function 
to a phase curve y, the total energy function can be converted into a function of 
q,p. One take H consisting of this new function and a constant as a Hamiltonian 
quantity, such that a Hamilton's system (i.e., a conservative system) is obtained. 
Under the initial condition mentioned above, the solution curve of the conserva- 
tive system is the same as that of the original nonconservative system; under other 
initial conditions the solution curve of the conservative is different from that of 
the original nonconservative system. Since one defines the forces JTl 12012 1 1221) in 
the new system, the Hamiltonian quantity of the conservative can be thought of as 
the mechanical energy of the new conservative system as Eq.i 
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The Hamiltonians of the new conservative systems in general are not analyti- 
cally integrable, unless the original mechanical system is integrable. The reason is 
that the work done by damping force depends on the phase curve. If the system is 
integrable, then the phase curve can be explicitly written out, the system has an an- 
alytical solution, and therefore the work done by damping force can be explicitly 
integrated. Subsequently, the Hamiltonian H can be explicitly expressed. Most 
systems do not have an analytical solution. Despite this, the Hamilton quantity, 
coordinates and momentum must satisfy Eq.(|4]l under a certain initial condition. 
Why had KleinlSj written, "Physicists can make use of these theories only very 
little, an engineers nothing at all"? The answer: when one is seeking an analytical 
solution to a classical mechanics problem by utilizing Hamiltonian formalism, in 
fact one must inevitably convert the problem back to Newtonian formalism. This 
means that an explicit form of Hamiltonian quantity is not necessary for classical 
mechanics. What is important is the relationship between q, p and the Hamiltonian 
quantity embodied in the Hamilton's Equation. 

According to the conclusion above, we can consider the Euler time-centered 
difference scheme for dissipative systems from a Hamiltonian perspective. 



3 Discussion on the Euler Time-Centered Difference Scheme 

3.1 Introducing the Euler Time-Centered Scheme for Conservative Systems 

Feng ll3] fT4ll 1 5 I l4l constructed a series of symplectic difference schemes via two 
approaches: the first approach discretizes Hamilton's equation and utilizes the 
property of the Cayley transform to demonstrate that the map ^ : z" — ?> z""*"' is 
a symplectic map; the second approach is a so-called generating function method. 
Feng had represented the first approach as below: 

Suppose // is a differentiable function of 2« variables p\,- ■ ■ , /?« , ^i , ■ ■ ' i ^«> 
the Hamilton's equations are represented as: 

p = -Hg, q=Hp, (25) 

where p = [pi,- ■ ■ ,pn],q = [qi,- ■ ■ ,qn], Hq = dH/dq and Hp = dH/dp. Let 
Z = [p,q] , and Eq. ( 125b can be further represented as: 

z = J-^H,, (26) 

where H, = [Hq,Hpf, 



J = 



O In 
-InO 



(27) 



where /„ and / denote unity matrices. If // is a quadratic form: 

H{z) = \z^Cz, C^ = C, (28) 

then canonical equations (l25 126b can be rewritten as: 

Bz, B = J-^C. (29) 



dz „ „ ,_i, 



df 

Paper ifTSJI gives the following definition 
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Definition 1 A matrix B of order 2n is called infinitesimal symplectic if 

jb + b'^j = o, 

where is a null matrix. All infinitesimal symplectic matrices form a Lie alge- 
bra sp{2n) with commutation operation [A, B] = AB ~ BA, and sp{2n) is the Lie 
algebra of Lie group Sp{2n) known as the symplectic group. 

According to this definition, 6 = J^'C in Eg. (129b can be considered as an 
infinitesimal symplectic matrix. In the paper|l3] a number of symplectic schemes 
for Eq.(|29ll were proposed, one of which is the Euler time-centered scheme: 

z"+^-z" „z"+^+z" 



B . (30) 

X 2 

The transition z" -^ z""*"' is given by the following linear transformation Fj which 
coincides with its own Jacobian 

F, = (/)(-|e) = (/-|e)-'(/ + |e). (3i) 

Paper ifTSJI gives the following theorem: 

Theorem 2 If B E sp{2n) and \l + B\ ^ 0, then F = {I + B)-\l ~ B) eSp{2n), 
the Cayley transform of B. 

According to Theorem |2l F-c can be considered as the Cayley transform of 
the infinitesimal symplectic matrix — t6/2, and consequently F-i is a symplectic 
matrix, and z" — > z"^^ is symplectic. 



3.2 Apply the Time-centered Difference Scheme Indirectly to Damping 
Oscillators 

Based on the discussion in Section |Z21 we have created a flow chart (Fig|2]) that 
shows the process of applying the time-centered difference scheme indirectly to 
the damping oscillatordlTll. We have defined a conservative system through an 
equivalent stiffness matrix represented by Eq.( l20b(l21l)(l22b . If one needs to formu- 
late the equivalent stiffness matrix in numerical schemes, one must have a numer- 
ical solution at time t^'^^, such that an equivalent elastic restoring force can be 
thought of as a function of (^^ +q^ )/2. Hence we first apply the time-centered 
difference scheme directly to the original dissipative system. Consequently the 
algorithm can be depicted by Fig. |2l 

The first step is to apply the time-centered difference scheme directly to Eq. (ll7b 

T ^ 2 
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Apply Time-centered Scheme to Obtain 
Probe Soiution 



Utilize Probe Solution To Compute 
Equivalent Stiffness 



Apply Time-centered Scheme to 
Corresponding Conservative System 



Fig. 2 Flowchart of applying the time-centered difference scheme indirectly to damping Oscil- 
lators in a time-step 



If we let d[(^'^ + ^*+^)/2]/df = (^'■'+' - ^'■')/t according to the definition of the 
time-centered difference scheme, Eq. (l32b can be represented as: 



„k+\ 



= F 



(33) 



where the transition matrix is 



[ ' r'] 


i 


.-<-' 





C) / 



(34) 



Obviously the difference scheme is not a symplectic scheme, because the transi- 
tion matrix F^ is not a symplectic matrix. The prerequisite proposition llSII for the 
proof to this point is as the following: 

Proposition 2 Matrix S = M'^ N e Sp{2n), iffMJM''' = NJN^ 



The proof of this point is given as follows: 
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Proof 


Since 














% 




J 


' 1 T(i./< + i.C)" 




= 


t2 1 1 1 

-{-■K+--C) + l 

T^ 1 1 




; 


' I'] 


J 






= 






and according to the Proposition |2] 


F^- 


[.(i-K + i-C) 


■ l' 

1 




1 


_-,i.K-i.c; 


T 

/ 


/" 



F\ is not a symplectic matrix. D 

Utilizing the scheme described by Eq.( l33l) we can obtain the numerical so- 
lution [^''+\ ;?''+'] which is taken as a probe solution in lieu of the analytical 
solution, and then execute the second step in Fig. |2l Let 



K = 



Kn 



K„. 



(35) 






^^})/(T(^r'+^n), 



where K is the numerical approximation of the equivalent stiffness matrix([ 
Thus we would obtain the numerical approximation of the conservative system([ 

Then we consider to apply the time-centered scheme to the conservative svstem.( l23l l. 
Suppose the solution vector at the time f*^"*"^ is z*^+^ = ["^^^^p^"*"'] According to 



the time-centered difference scheme defined by Eg. ( 1281 
centered scheme for the conservative svstem( l23] l: 



, we have a time- 



~qk+i. 



pk+i_pk 



pk+l^pk 
2 

~iK+K) 



~A:+1 



(36) 
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The time-centered scheme above can be represented as 



pk+i 



where 



(37) 



/ 



:iK + K) 



I 



-■I 



-■I 



'-{K+K) I 



(38) 



According to Eq.dJTJ) and Ea.([38]). the map z'' 

Ff must be a symplectic matrix. The proof of this point is as follows 



z*^+' must be symplectic, because 



Proof 



' T 







i' 



:{K+K) 



-1 



Since 



/- 
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K) 
and according to Definition [T] 
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is an infinitesimal symplectic matrix. Therefore, according to Theorem |2] F^ is a 
symplectic matrix. 



Finally, by substituting ^ = \q'^ , p*^] and Eqs.(l35Tl into Eg.dJTI) and Eq.(l38]). we 
can obtain the numerical solution z*^+^. Then if we repeat the process that consists 
of Eg. ( 1331) ( |35]) ( 137b . we would get a series of numerical solutions for Eg. (117b that 
are canonical for its substituting conservative system. 
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Fig. 3 Total energy behavior of integrators for an one-dimensional damping oscillator 



From the derivation above, we can obtain the numerical solution z^ by direct 
application to the dissipative svstem( ll7b and the numerical solution l*""^' by indi- 
rect application. Although F^ ^ F^, from Eq.([33l), Eq.®, Eq.lO and Eq.(l36l) 
we can derive 



f,k+\ 



jk+l 
n<:+l 



V - ^-Kq^ + IxCq^ + 4Tp^ 



4 + t2/< 
AxKq'^ + x-Kp^ 



2tC 
V2TCp''-Ap^ 



4 + t2K' + 2tC 



(39) 



According to Eq.(|39]l, the result of the transition z*^ -^ 2*^+' is exactly identical 

It is a known fact, that the 



.k+\ 



k+\ 



to that of z* — > z*""^', and hence we have z' 

linear transition matrix between vectors z'' and z"+' should not be unique; hence 
F^ may not be equal to F^. Therefore, we can say that the time centered differ- 
ence scheme(l33b for a dissipative mechanical system is in fact also a symplectic 
scheme for the substituting conservative system. Consequently, one does not need 
to execute the second and third steps in Fig. |2] 



3.3 Numerical Examples 

Although the time-centered difference scheme is used widely, the total energy be- 
havior of this scheme for dissipative systems has not been illustrated via numerical 
examples. We now have two numerical examples to show the total energy behavior 
for damping systems. 

The first numerical example involves an one-dimensional damping oscillator 



-2x + 0.05i- = 
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Fig. 4 Total energy behavior of integrators for a two-dimensional damping oscillator 
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Fig. 5 1st displacement of a two-dimensional damping oscillator 



with initial values jcq = 0.1, io = 0.2. The time-centered scheme with step- 
size=0.2 and the 4th-Runge-Kutta method with step-size=0.2 are used to compute 
the one-dimensional problem. The total energy result is shown in Fig. [3] 

The second numerical example involves a two-dimensional damping oscillator 

x\ + 3x1 + 0.03x1 - 0.01x2 = 
X2 + 3X2 - O.Olx'i + 0.01X2 = 0, 

with initial values xi|;=o = 0.1, xi|(=o=0.1, X2|r=o=0.2, X2|r=o = 0.2. The 
same method as earlier has been employed to compute this problem. The numer- 
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ical result of the 1st displacement is shown in Fig|5] and the total energy result in 
Fig. SI 

As shown in Fig. [3] and Fig.|4l the Euler time-centered difference scheme can 
preserve the total energy, the primary reason being is that the foundation of the 
Euler time-centered difference scheme is a Hamilton's system, the Hamiltonian 
quantity of which is the sum of the total energy and a constant. By comparing Fig. 
|4]and Fig.|5] one will find that the period of the result of the time-centered scheme 
is shorter than that of the result of the Runge-Kutta method. The reason is clear: the 
Runge-Kutta method would cause artificial energy growth. This result can be con- 
sidered as a generalization of the conclusion in the paper lll6l . which asserts that 
the Euler time-centered scheme preserves the mechanical energy of conservative 
oscillators. By comparing the results in this paper and the work in the paperQ, we 
can identify the difference. As presented in Fig. 6.1 and Fig. 6.2 in the paper Q, 
'The variational integrators simulate energy decay, unlike standard methods such 
as Runge-Kutta.', there is no accurate criterion of numerical integration methods 
for dissipative system just as the accurate criterion for conservative systems which 
is a horizontal line(see Fig. 4.1 and Fig. 4.2 in the paperQ), its rationale has have 
been presented in Sec(T] 

4 conclusions 

We can conclude that a dissipative mechanical system has such properties: for any 
nonconservative classical mechanical system and any initial condition, there ex- 
ists a conservative one, the two systems share one and only one common phase 
curve; the Hamiltonian of the conservative system is the sum of the total energy 
of the nonconservative system on the aforementioned phase curve and a constant 
depending on the initial condition. We can further conclude, that a dissipative 
problem can be reformulated as an infinite number of non-dissipative problems, 
one corresponding to each phase curve of the dissipative problem. One can avoid 
having to change the definition of the canonical momentum in the Hamilton for- 
malism, because under a certain initial condition the motion of one of the group 
of conservative systems is the same as the original dissipative system. 

From a Hamiltonian perspective, this paper has revealed that the transition 
matrix of z^ -^ Z*^^' is not unique, the transition matrix may be non-symplectic or 
symplectic, the numerical solution of the original dissipative obtained through the 
time-centered difference scheme system is equal to that of the substituting con- 
servative system obtained through a time-centered difference scheme. In addition, 
we have discovered that the time-centered scheme preserves the total energy of 
dissipative systems. Based on the above, we might explain why the time-centered 
difference scheme is more accurate than an unsymplectic one for dissipative prob- 
lem. For this reason we might be able to choose and construct better algorithms. 
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